!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

module mMixture

use mTypes

implicit none

type tTie_Line
  real(rb), allocatable :: Z(:), X(:), Y(:)
  contains
    procedure :: read => Tie_Line_read
end type tTie_Line

type tCondition
  real(rb) :: Temperature
  integer  :: Number_of_Tie_Lines
  type(tTie_Line), allocatable :: Tie_Line(:)
  character(sl) :: Reference
  contains
    procedure :: Read => tCondition_Read
end type tCondition

type tMixture
  integer :: Number_of_Components
  character(sl), allocatable :: Component(:)
  integer :: Number_of_Conditions
  type(tCondition), allocatable :: Condition(:)
  contains
    procedure :: Read => tMixture_Read
end type tMixture

contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine Tie_Line_read( tie_line, unit, NC )
    class(tTie_Line), intent(inout) :: tie_line
    integer,          intent(in)    :: unit, NC
    if (allocated(tie_line % Z)) deallocate( tie_line % Z, tie_line % X, tie_line % Y )
    allocate( tie_line % Z(NC), tie_line % X(NC), tie_line % Y(NC) )
    read(unit,*) tie_line % Z, tie_line % X, tie_line % Y
  end subroutine Tie_Line_read

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tCondition_Read( condition, unit, NC )
    class(tCondition), intent(inout) :: condition
    integer,           intent(in)    :: unit, NC
    integer :: i
    read(unit,*); read(unit,*) condition % Temperature
    read(unit,*); read(unit,*) condition % Number_of_Tie_Lines
    allocate( condition % Tie_Line(condition % Number_of_Tie_Lines) )
    read(unit,*)
    do i = 1, condition % Number_of_Tie_Lines
      call condition % Tie_Line(i) % read( unit, NC )
    end do
    read(unit,*); read(unit,*) condition % Reference
  end subroutine tCondition_Read

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tMixture_Read( mixture, unit )
    class(tMixture), intent(inout) :: mixture
    integer,         intent(in)    :: unit
    integer :: i
    read(unit,*); read(unit,*) mixture % Number_of_Components
    allocate( mixture % Component(mixture % Number_of_Components) )
    read(unit,*)
    do i = 1, mixture % Number_of_Components
      read(unit,*) mixture % Component(i)
    end do
    read(unit,*); read(unit,*) mixture % Number_of_Conditions
    allocate( mixture % Condition( mixture % Number_of_Conditions) )
    do i = 1, mixture % Number_of_Conditions
      call mixture % Condition(i) % Read( unit, mixture%Number_of_Components )
    end do
  end subroutine tMixture_Read

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

end module mMixture

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

module mModel

use mTypes

implicit none

type tSubgroup
  integer      :: Index
  character(8) :: Name
  integer      :: Main
  character(8) :: Main_Name
  real(rb)     :: R
  real(rb)     :: Q
  integer,  allocatable :: Incidence(:)
end type tSubgroup

type, abstract :: tModel
  integer :: Number_of_Components
  integer :: Number_of_Subgroups
  integer :: Number_of_Conditions
  type(tSubgroup), allocatable :: Subgroup(:)

  real(rb), allocatable :: R(:), Q(:), Omega(:,:)

  contains
    procedure :: Identify_Subgroups => tModel_Identify_Subgroups
    procedure :: Read_Group_Parameters => tModel_Read_Group_Parameters
    procedure :: Read_Interaction_Parameters => tModel_Read_Interaction_Parameters
    procedure(tModel_Allocate_Parameters), deferred :: Allocate_Parameters
    procedure(tModel_Get_Parameters), deferred :: Get_Parameters
    procedure :: Define_Matrices => tModel_Define_Matrices
end type tModel

abstract interface

  subroutine tModel_Allocate_Parameters( this )
    import :: tModel
    class(tModel), intent(inout) :: this
  end subroutine tModel_Allocate_Parameters

  subroutine tModel_Get_Parameters( this, unit, i, j )
    import :: tModel
    class(tModel), intent(inout) :: this
    integer,       intent(in)    :: unit, i, j
  end subroutine tModel_Get_Parameters

end interface

contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tModel_Identify_Subgroups( this, component, model )
    class(tModel), intent(inout) :: this
    character(*),  intent(in)    :: component(:), model
    integer :: i, N, j, Nrep, Indx
    logical :: Found
    type Item
      integer :: Index
      integer, allocatable :: Nrep(:)
      type(Item), pointer :: Next
    end type Item
    type(Item), pointer :: First => null(), Current
    associate ( NC => this % Number_of_Components, NSG => this % Number_of_Subgroups )
      NC = size(component)
      NSG = 0
      do i = 1, NC
        open( unit = 10, file = "compounds/"//trim(component(i))//"."//trim(model), status = "old" )
        read(10,'(/,I10,/)') N
        do j = 1, N
          read(10,*) Nrep, Indx
          Current => First
          Found = .false.
          do while (associated(Current).and.(.not.Found))
            Found = Current % Index == Indx
            if (Found) Current % Nrep(i) = Nrep
            Current => Current % Next
          end do
          if (.not.Found) then
            allocate( Current )
            Current % Index = Indx
            allocate( Current % Nrep(NC) )
            Current % Nrep = 0
            Current % Nrep(i) = Nrep
            Current % Next => First
            First => Current
            NSG = NSG + 1
          end if
        end do
        close( 10 )
      end do
      allocate( this % Subgroup( NSG ) )
      Current => First
      j = NSG
      do while (associated(Current))
        this % Subgroup(j) % Index = Current % Index
        allocate( this % Subgroup(j) % Incidence(NC) )
        this % Subgroup(j) % Incidence(:) = Current % Nrep
        Current => Current % Next
        j = j - 1
      end do
    end associate
  end subroutine tModel_Identify_Subgroups

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tModel_Read_Group_Parameters( this, model )
    class(tModel), intent(inout) :: this
    character(*),  intent(in)    :: model
    integer :: i, error
    logical :: Found
    type(tSubgroup) :: SG
    open( unit = 10, file = "parameters/"//trim(model)//"_groups.dat", status = "old" )
    do i = 1, this % Number_of_Subgroups
      rewind(10)
      read(10,*)
      read(10,*)
      Found = .false.
      error = 0
      do while (.not.Found.and.(error == 0))
        read(10,'(I5,A8,I14,A8,F19.0,F12.0)',iostat=error) &
          SG % Index, SG % Name, SG % Main, SG % Main_Name, SG % R, SG % Q
        if (error == 0) Found = SG % Index == this % Subgroup(i) % Index
      end do
      if (Found) then
        this % Subgroup(i) % Name = SG % Name
        this % Subgroup(i) % Main = SG % Main
        this % Subgroup(i) % Main_Name = SG % Main_Name
        this % Subgroup(i) % R = SG % R
        this % Subgroup(i) % Q = SG % Q
      else
        write(*,'("Error: subgroup ",I3," not found in file ",A)') &
          i, "parameters/"//trim(model)//"_groups.dat"
        stop
      end if
    end do
    close( 10 )
  end subroutine tModel_Read_Group_Parameters

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tModel_Read_Interaction_Parameters( this, model )
    class(tModel), intent(inout) :: this
    character(*),  intent(in)    :: model
    integer  :: i, j, imain, jmain, first, second, ii, jj, io
    open( unit = 10, file = "parameters/"//trim(model)//"_parameters.dat", status = "old")
    call this % Allocate_Parameters()
    do i = 1, this % Number_of_Subgroups
      imain = this % Subgroup(i) % Main
      do j = i + 1, this % Number_of_Subgroups
        jmain = this % Subgroup(j) % Main
        if (jmain /= imain) then
          first = min(imain,jmain)
          second = max(imain,jmain)
          rewind( 10 )
          read(10,'(/)')
          ii = 0
          jj = 0
          io = 0
          do while ((io == 0).and.( (ii < first).or.( (ii == first).and.(jj < second) ) ))
            read(10,*)
            read(10,'(I7,I8)',advance="no",iostat=io) ii, jj
          end do
          if ((ii == first).and.(jj == second)) then
            if (imain < jmain) then
              call this % Get_Parameters( 10, i, j )
            else
              call this % Get_Parameters( 10, j, i )
            end if
          else
            write(*,'("WARNING: interaction parameters between subgroups ",I3,'// &
                    '" and ",I3," were not found in file",A,".")') & 
                    i, j, "parameters/"//trim(model)//"_parameters.dat"
          end if
        end if
      end do
    end do
    close(10)
  end subroutine tModel_Read_Interaction_Parameters

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tModel_Define_Matrices( this )
    class(tModel), intent(inout) :: this
    integer :: NC, NSG, i, j
    real(rb) :: Rg(this%Number_of_Subgroups), Qg(this%Number_of_Subgroups)
    real(rb) :: Nu(this%Number_of_Components,this%Number_of_Subgroups)
    NC = this%Number_of_Components
    NSG = this%Number_of_Subgroups
    Rg = [(this%Subgroup(i)%R, i = 1, NSG)]
    Qg = [(this%Subgroup(i)%Q, i = 1, NSG)]
    forall (i = 1:NC, j = 1:NSG)
      Nu(i,j) = this % Subgroup(j) % Incidence(i)
    end forall
    if (allocated(this%R)) deallocate( this%R, this%Q, this%Omega )
    allocate( this%R(NC), this%Q(NC), this%Omega(NC,NSG) )
    this%R = matmul( Nu, Rg )
    this%Q = matmul( Nu, Qg )
    forall (i = 1:NC) this%Omega(i,:) = Nu(i,:)*Qg
  end subroutine tModel_Define_Matrices

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

end module mModel

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

module mUNIFAC

use mModel

implicit none

type, extends(tModel) :: tUNIFAC
  real(rb), allocatable :: A(:,:)
  contains
    procedure :: Allocate_Parameters => tUNIFAC_Allocate_Parameters
    procedure :: Get_Parameters => tUNIFAC_Get_Parameters
end type tUNIFAC

contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tUNIFAC_Allocate_Parameters( this )
    class(tUNIFAC), intent(inout) :: this
    integer :: NSG
    NSG = this % Number_of_Subgroups
    if (allocated(this % A)) deallocate( this % A )
    allocate( this % A(NSG,NSG) )
    this % A = 0.0_rb
  end subroutine tUNIFAC_Allocate_Parameters

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  subroutine tUNIFAC_Get_Parameters( this, unit, i, j )
    class(tUNIFAC), intent(inout) :: this
    integer,        intent(in)    :: unit, i, j
    read(10,*) this % A(i,j), this % A(j,i)
  end subroutine tUNIFAC_Get_Parameters

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

end module mUNIFAC

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

program UNIFAC

use mMixture
use mUNIFAC
use Fortrix

implicit none

integer :: i, j

type(tMixture) :: Mixture
type(tUNIFAC) :: Model

type(scalar) :: T
type(matrix) :: r, q, m, Omega, A, vec_1, q4
type(matrix) :: G, eps
type(matrix) :: x, phi, theta, L
type(matrix) :: ln_gamma

call Mixture % Read( 5 )
call Model % Identify_Subgroups( Mixture % Component, "unifac" )
call Model % Read_Group_Parameters( "unifac" )
call Model % Read_Interaction_Parameters( "unifac_lle" )

write(*,'("No.  Subgroup      Main       R       Q")')
do i = 1, Model % Number_of_Subgroups
  associate( SG => Model % Subgroup(i) )
    write(*,'(I3,2A10,2F8.4)') i, trim(SG%Name), trim(SG%Main_Name), SG%R, SG%Q
  end associate
end do

call fortrix_startup()
call Model % Define_Matrices()

r = new_matrix( Model % R )
q = new_matrix( Model % Q )
q4 = new_scalar(4.d0)*q
m = ones( Model % Number_of_Components ) - (q4 + q)
Omega = new_matrix( Model % Omega )
A = new_matrix( Model % A )
vec_1 = ones( Model % Number_of_Subgroups )

do i = 1, Mixture % Number_of_Conditions
  T = new_scalar( Mixture % Condition(i) % Temperature )
  G = Omega*exp(-inv(T)*A)
  eps = q + (Omega.o.G)*vec_1
  do j = 1, Mixture % Condition(i) % Number_of_Tie_Lines
    x = new_matrix( Mixture % Condition(i) % Tie_Line(j) % X )
    phi = r*inv(.t.r * x)
    theta = q*inv(.t.q * x)
    L = G*Dinv(.t.G * x)
    ln_gamma = D(m)*log(phi) - D(q4)*log(theta) - phi*(.t.m)*x + m + eps - Omega*log(.t.G * x) - L*(.t.Omega)*x
    print*
    call print_matrix( ln_gamma, "F11.4" )
  end do
end do

call fortrix_shutdown()

end program UNIFAC
